% !Mode:: "TeX:UTF-8" 

\BiChapter{基于概率拓扑的构型规划方法}{Configuration planning under belief space with probabilistic network topology}

\BiSection{引言}{Introduction}

在连续信念空间的构型规划框架中，多智能体系统需要对每一个候选控制输入计算它们所对应的目标函数，然后对比评估该控制输入的表现优劣。
因此式~\eqref{eq. ACTCO_objective_function}的计算需要在当前规划时刻$k$获取规划期限$k+1:k+L$内的观测数据。
然而，在规划任务中，对未来状态的预测并不存在实际的测量操作，因此规划期限内的观测是未知的。
前一章采用了最大似然假设，即使用最大似然观测来近似表征未来规划期限内的观测数据。然而这种方法仅评估了最大似然这一种未来分支，并未评估所有可能性。

当考虑观测设备存在的实际约束时，观测连通变量的具体分布由多智能体节点的预测状态以及装备传感器的观测模型所决定，因此实际应用中多智能体网络中的观测拓扑将保持时变。由于多智能体系统受到来自运动过程和观测过程的噪声干扰，因此多智能体系统在未来规划期限内的预测状态是随机分布的，从而导致未来时刻的观测数据在规划时刻也是随机分布的。

已有研究成果的一些定性分析表明，多智能体内部的观测数量越多，协同定位的精度也会越高。因此本章将聚焦于预测多智能体系统在规划期限内的观测拓扑的变化。与实际观测变量的连续分布所不同的是，观测连通变量是一类伯努利变量，只有$0$或$1$两种状态，因此多智能体系统在规划期限内所有可能的观测拓扑分支的数量是有限的。因此，如果可以根据规划时刻已知信息，准确预测规划周期内的观测拓扑概率，则可以提升构型规划方法对系统未来行为预测准确性，进一步提升构型规划方案对环境不确定性的鲁棒性。

目前，已有的一些信念空间规划工作\citeup{platt2010belief,indelman2015planning,pathak2018unified}中对拓扑概率的计算主要是三类近似方法，即采用最大似然观测假设的伯努利模型，线性近似模型和抽样频率近似模型。这些方法全都是经验近似的概率计算方法，计算精度没有理论保障。另外，目前的研究成果缺乏针对利用精确的拓扑概率提升信念空间规划性能的深入研究。

针对以上问题，本章研究一类拓扑概率的精确计算方法，并以此为基础设计了一类基于概率拓扑的构型规划方法。本章不再仅使用最大似然假设近似计算观测连通变量，而是在观测模型的基础上，通过对节点之间相对观测过程的精确建模计算观测连通变量的分布。本章在高斯噪声模型假设下，针对一类圆盘通信测量模型，采用正态二次型分布来建模两个节点在未来时刻的观测连通变量。基于正态二次型多项式展开定理，提出一类计算精度有理论保障的观测连通概率自适应计算方法；然后将该算法应用到信念空间构型规划框架，提出基于概率拓扑的目标函数评估方法，以提升了构型规划方法对多智能体系统未来行为预测的准确性。最后，通过数值仿真和对比验证了拓扑概率计算方法的准确性，以及基于概率拓扑的构型规划方法的有效性。


%本章研究一类基于路径规划的主动协同定位方法。已有的研究成果表明，当系统可以获取绝对定位信息时（例如系统中至少有一个节点存在GNSS测量信息，或者可以观测到一个固定的锚节点），则该系统的协同定位过程是可观测的，因此系统中所有节点的稳态定位误差是有界的。该界限与节点初始的位置误差无关，而是受到系统内部的相对测量拓扑和节点装配的对内/对外测量传感器的精度的影响。因此作为一个可控变量，系统每一时刻的相对观测拓扑可以被用来优化协同定位的性能。另外，当系统内每个智能体装备的测量和通信设备的最大观测和通信距离受限时，系统内部的通信观测拓扑将不再是固定不变的，而是会随着节点之间相对位置的变化而改变。又由于各个智能体都存在自定位误差，因此主动协同定位过程对未来时刻系统内部通信测量拓扑的预测将是随机的、不确定的。
%
%基于以上问题，本章首先基于优化理论，将每个智能体的运动轨迹视作优化变量，控制多智能体内部相互观测构型，以获取最佳定位数据，通过求取最优的运动路径，提升协同定位任务的性能。其次，本章基于前一章的研究基础，引入概率拓扑的概念，即将系统在未来时刻的通信和观测拓扑视作随机和不确定的，从而设计了一类基于概率拓扑的主动协同定位方法。该方法通过前一章设计的未来连通概率预测算法，计算未来时刻每一种可能的观测拓扑出现的概率，从而预测协同定位任务在给定路径上性能表现的期望值。以此来评价每条路径的优先值，从而规划最佳的路径输入。


\BiSection{概率拓扑对目标函数的影响}{The impact of probabilitic topology on objective function}

本节将对式~\eqref{eq. ACTCO_objective_function}目标函数中的取期望操作进行展开。根据取期望操作的定义将式~\eqref{eq. ACTCO_objective_function}转变为积分操作，从而分解目标函数的计算过程，并重点分析概率拓扑对目标函数计算的影响。首先，规划期限内的观测变量也包含实际观测数据和观测拓扑两部分，即$\vect{Z}^{k+1:k+L}=\{ \vect{\vec{Z}}^{k+1:k+L},\mathcal{E}^{k+1:k+L} \}$，因此可以将式~\eqref{eq. ACTCO_objective_function}中的取期望操作具化为积分操作，并使用全概率公式得：
\begin{equation*} 
	\begin{aligned}
		& J_{k}(\vect{U}^{k:k+L-1}) =   \int_{\vect{\vec{Z}}^{k+1:k+L}} \int_{\mathcal{E}^{k+1:k+L}} p( \vect{\vec{Z}}^{k+1:k+L}, \mathcal{E}^{k+1:k+L} |\mathcal{H}^{0:k+L|k} ) \\
		& \left[ \sum_{l=1}^{L} c^{l} \left( p(\vect{X}^{0:k+l}|\mathcal{H}^{0:k+l|k}, \vect{\vec{Z}}^{k+1:k+l}, \mathcal{E}^{k+1:k+l}), \vect{U}^{k:k+l-1} \right)  \right]  \\
		& = \int_{\mathcal{E}^{k+1:k+L}} p( \mathcal{E}^{k+1:k+L} |\mathcal{H}^{0:k+L|k} )  \int_{\vect{\vec{Z}}^{k+1:k+L}}  p( \vect{\vec{Z}}^{k+1:k+L}  |\mathcal{H}^{0:k+L|k},\mathcal{E}^{k+1:k+L} ) \\
		& \left[ \sum_{l=1}^{L} c^{l} \left( p(\vect{X}^{0:k+l}|\mathcal{H}^{0:k+l|k}, \vect{\vec{Z}}^{k+1:k+l}, \mathcal{E}^{k+1:k+l}), \vect{U}^{k:k+l-1} \right)  \right] 
	\end{aligned}	
\end{equation*}
由于在规划期限内，观测、通信拓扑的分布是离散的，其值域$\vect{E}$的规模是有限的。给定集群的规模为$N$，环境中存在$M$个锚节点，规划期限为$L$，则未来规划期限内拓扑集合$\vect{E}$中最多可以包含$N_{con}$个元素，
$$N_{con} = 2^{N_s}, N_s = \left[ \frac{L(M+N)(M+N-1)}{2} \right]$$

目标函数可进一步转变为：
\begin{equation} \label{eq. ACTCO_objective function1}
	\begin{aligned}
		& J_{k}(\vect{U}^{k:k+L-1}) =  \sum_{\mathcal{E}^{k+1:k+L}_{iter} \in \vect{E}} \underbrace{ p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k} ) }_{A} \int_{\vect{\vec{Z}}^{k+1:k+L}}  \underbrace{ p( \vect{\vec{Z}}^{k+1:k+L}  |\mathcal{H}^{0:k+L|k},\mathcal{E}^{k+1:k+L}_{iter} ) }_{B} \\
		& \left[ \sum_{l=1}^{L} c^{l} \left( \underbrace{ p(\vect{X}^{0:k+l}|\mathcal{H}^{0:k+l|k}, \vect{\vec{Z}}^{k+1:k+l}, \mathcal{E}^{k+1:k+l}_{iter})}_{C}, \vect{U}^{k:k+l-1} \right)  \right] 
	\end{aligned}	
\end{equation}
因此，当给定一组待评估的候选控制${\vect{U}}^{k:k+L-1}$，其所对应的目标函数包含三部分：$A$为未来时刻不确定性拓扑的似然项；$B$为在给定观测拓扑前提下，观测传感器实际测量数据的似然；$C$为在给定观测拓扑和实际观测数据前提下，对联合状态向量后验概率密度分布的预测。下面分别对$A,B,C$三项进行详细分析：

(1) \textbf{B-测量似然}$p( \vect{\vec{Z}}^{k+1:k+L}  |\mathcal{H}^{0:k+L|k},\mathcal{E}^{k+1:k+L}_{iter} )$:

首先分析测量似然项，由于该项并非本章所研究的重点，为了简化问题并突出概率拓扑的影响，本文将采用最大似然假设来计算测量似然项的取值，即
$$p( \vect{\vec{Z}}^{k+1:k+L}  |\mathcal{H}^{0:k+L|k},\mathcal{E}^{k+1:k+L}_{iter} ) = \left\lbrace \begin{array}{ll}
	1, & \quad \text{如果 } \vect{\vec{Z}}^{k+1:k+L} = \vect{\vec{Z}}^{k+1:k+L}_{ML} \\
	0, & \quad \text{否则}
\end{array} \right.
$$
因此，目标函数~\eqref{eq. ACTCO_objective function1}可转换为：
\begin{equation} \label{eq. ACTCO_objective function2}
	\begin{aligned}
		& J_{k}(\vect{U}^{k:k+L-1}) =  \sum_{\mathcal{E}^{k+1:k+L}_{iter} \in \vect{E}} \underbrace{ p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k} ) }_{A} \\
		& \left[ \sum_{l=1}^{L} c^{l} \left( \underbrace{ p(\vect{X}^{0:k+l}|\mathcal{H}^{0:k+l|k}, \vect{\vec{Z}}^{k+1:k+l}_{ML}, \mathcal{E}^{k+1:k+l}_{iter})}_{C}, \vect{U}^{k:k+l-1} \right)  \right] 
	\end{aligned}	
\end{equation}

依据最大似然假设，实际观测数据$\vect{\vec{Z}}^{k+1:k+L}_{ML}$可由式~\eqref{eq. ACOTO_priorBeliefGivenControl}中的先验信念分布生成。因此，在后文的具体算例中，需要首先调用附录~\ref{ch5sec. appendix}中的估计引擎，求解符合式~\eqref{eq. ACOTO_priorBeliefGivenControl}的高斯分布状态均值$\vect{\bar{X}}^{k+1:k+L|k}$，然后使用该状态均值和测量设备的观测模型生成无误差的最大似然观测数据$\vect{\vec{Z}}^{k+1:k+L}_{ML}$，每一组数据的存在与否受给定拓扑$\mathcal{E}^{k+1:k+L}_{iter}$的约束。

(2) \textbf{C-未来状态的后验概率密度分布}$p(\vect{X}^{0:k+l}|\mathcal{H}^{0:k+l|k}, \vect{\vec{Z}}^{k+1:k+l}_{ML},\mathcal{E}^{k+1:k+l}_{iter})$：

该项表示在假定的通信观测拓扑$\mathcal{E}^{k+1:k+l}_{iter}$以及实际观测数据$\vect{\vec{Z}}^{k+1:k+l}_{ML}$下，系统在未来规划期限内状态向量的后验概率密度分布。由于本文采用高斯噪声假设，因此计算该后验概率密度分布等价于寻找一组高斯分布的状态均值$\vect{\bar{X}}^{0:k+l}$和分布协方差矩阵$\Sigma_{\bar{X}^{0:k+l}}$。
因此该项的计算与估计问题相似，可直接调用附录~\ref{ch5sec. appendix}中的滤波引擎和光滑引擎进行求解。

(3) \textbf{A-拓扑似然}$p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k} )$:

针对拓扑似然概率项$p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k} )$，在系统在未来时刻的联合状态变量$\vect{X}^{k+1:k+L}$上再次运用全概率公式，式~\eqref{eq. ACTCO_objective function1}中的$A项$可作如下展开：
\begin{equation*}
	\begin{aligned}
		p( \mathcal{E}^{k+1:k+L}_{iter} &|\mathcal{H}^{0:k+L|k} )  = \int_{\vect{X}^{k+1:k+L}} p( \mathcal{E}^{k+1:k+L}_{iter}, \vect{X}^{k+1:k+L} |\mathcal{H}^{0:k+L|k} ) \\
		& = \int_{\vect{X}^{k+1:k+L}} p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k},  \vect{X}^{k+1:k+L} ) p(\vect{X}^{k+1:k+L} | \mathcal{H}^{0:k+L|k} )
	\end{aligned}
\end{equation*}

本章假设未来网络拓扑集合$\mathcal{E}^{k+1:k+L}$中的任意两个元素在空间上和时间上都是互相独立的。传统的状态估计问题通常直接使用马尔可夫独立假设，即确定了系统的当前状态后，系统在未来时刻的行为与系统的过去状态无关。该假设认为当前时刻系统的状态是完备的，即过去时刻的任何信息，包括系统状态、控制输入以及观测数据等，对未来的系统行为不会产生影响。因此，马尔可夫独立假设保证了测量数据在时间上的相互独立性。另外，在某一时刻，各个节点在执行测量操作时，通常认为每个节点的测量操作是独立进行的，因此保证了测量数据在空间上的独立性。由此可以推断，在状态估计问题中，网络拓扑集合$\mathcal{E}^{0:k}$中的任意两个元素在时间上和空间上都是互相独立的。本章主要考虑规划问题。由于在规划时刻并不存在未来的观测数据，本章将沿用估计问题中拓扑独立的性质，假设未来时刻的网络拓扑集合$\mathcal{E}^{k+1:k+L}$中任意两个元素在空间上和时间上也都是互相独立的，即
\begin{itemize}
	\item[(1)] 同一个时刻中的任意两条观测和通信连通的存在关是相互独立的，即
	\begin{equation} \label{eq. connectionSpatialIndependence}
		p(\mathcal{E}^{t}) = \prod_{i \in \mathcal{V}_R} \prod_{j \in \vect{N_i^k}} p(s_{i,j}^t), \forall  k+1 \le t \le k+L
	\end{equation}
	\item[(2)] 任意两个时刻之间系统的拓扑关系是相互独立的，即
	\begin{equation} \label{eq. connectionTemporalIndependence}
		p(\mathcal{E}^{k+1:k+L}) = \prod_{t=k+1}^{k+L} p(\mathcal{E}^{t})
	\end{equation}	
\end{itemize}
一般而言，在保证同一时刻各个测量操作之间不存在相互干扰的情况下，规划过程测量数据的空间独立性将依然成立。然而由于估计问题中测量数据之间在时间上的独立性是基于马尔可夫假设推导而出的，因此时间独立性在规划过程的成立与否将取决于马尔科夫假是否依然有效。规划期限内马尔可夫属性的合理性有待进一步研究。

由拓扑概率的时空独立假设，目标函数\eqref{eq. ACTCO_objective function1}中与网络拓扑相关的概率项$p( \mathcal{E}^{k+1:k+L} |\mathcal{H}^{0:k+L|k} )$可使用全概率公式进一步展开为：

\begin{equation}\label{eq. PROB_decompositionConnPriori}
	\begin{aligned}
		p( \mathcal{E}^{k+1:k+L}_{iter} &|\mathcal{H}^{0:k+L|k} )  = \int_{\vect{X}^{k+1:k+L}} p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k},  \vect{X}^{k+1:k+L} ) p(\vect{X}^{k+1:k+L} | \mathcal{H}^{0:k+L|k} ) \\
		& = \int_{\vect{X}^{k+1:k+L}} \left[ \prod_{s_{ij}^t \in \mathcal{E}^{k+1:k+L}} p(s_{i,j}^t | \vect{p}_i^t, \vect{X}_j^t)  \right] p(\vect{X}^{k+1:k+L} | \mathcal{H}^{0:k+L|k} )
	\end{aligned}
\end{equation}
上式对多智能体节点未来状态向量$\vect{X}^{k+1:k+L}$的积分意味着，任意两个节点在未来时刻观测连通的概率是与它们节点位置的全分布$b(\vect{p}_i^t),b(\vect{X}_j^t)$相关的。即要确定$p(s_{i,j}^t | \vect{p}_i^t, \vect{X}_j^t)$的值，必须从状态向量的未来信念分布$p(\vect{X}^{k+1:k+L} | \mathcal{H}^{0:k+L|k})$中提取完整的节点分布信息。以高斯分布为例，不仅需要预测节点位置在未来规划期限内的位置均值$\vect{\bar{p}}_i^t,b(\vect{\bar{X}}_j^t$，还需要获取节点位置分布的协方差矩阵$\vect{\Sigma}_{\vect{p}_i^t}$和$\vect{\Sigma}_{\vect{X}_j^t}$。由于上式中全概率积分操作可以针对状态向量$\vect{X}^{k+1:k+L}$的任意分布，因此本章采用式~\eqref{eq. ACOTO_priorBeliefGivenControl}的先验概率分布。由此，式~\eqref{eq. PROB_decompositionConnPriori}可进一步拓展为：
\begin{equation}\label{eq. PROB_decompositionConnPriori1}
	\begin{aligned}
		p( \mathcal{E}^{k+1:k+L}_{iter} &|\mathcal{H}^{0:k+L|k} )  = \int_{\vect{X}^{k+1:k+L}} \left[ \prod_{s_{ij}^t \in \mathcal{E}^{k+1:k+L}} p(s_{i,j}^t | \vect{p}_i^t, \vect{X}_j^t)  \right] b(\vect{X}^{k+1:k+L|k})
	\end{aligned}
\end{equation}

后文的主要任务是当给定历史数据$\mathcal{H}^{0:k+L|k}$时，依据式~\eqref{eq. PROB_decompositionConnPriori1} 计算给定拓扑的概率值$p( \mathcal{E}^{k+1:k+L}_{iter} |\mathcal{H}^{0:k+L|k} )$。其中节~\ref{subsec. ch5sec3}设计一类针对圆盘观测模型的连通概率$p(s_{i,j}^t | \vect{p}_i^t, \vect{X}_j^t)$的计算方法。最后本章将这些分析结果应用到目标函数~\eqref{eq. ACTCO_objective function1}的计算过程中，从而提升主动协同定位构型的规划效率。

\BiSection{基于自适应参数选择的有限项近似概率计算方法}{Finite-term approximation for probability using adaptive power series extension} \label{subsec. ch5sec3}

本节将对一类圆盘观测传感器的观测连通模型进行建模，并设计一类概率预测算法，在两个邻接节点的分布已知时，计算他们之间存在连通的概率。在给出具体分析内容前，首先给出后文将使用的一个关于正态二次型分布的重要定理\citeup{provost1992quadratic}：
\begin{theorem}  \label{th.PROB_PowerSeriesExpansion}
	给定一个$p$维的多元正态分布向量$\vect{X} \sim \mathcal{N}(\vect{\mu}, \vect{\Sigma})$, $\vect{\Sigma} > 0$，定义它的二次型$Y = \vect{X}^T\vect{A}\vect{X}$，其中矩阵$\vect{A}$式一个正定对称矩阵，即$\vect{A} = \vect{A}^T > 0$。用$\vect{\lambda} = [\lambda_1,...,\lambda_p]^T$表示矩阵$\vect{\Sigma}^{\frac{1}{2}}A\vect{\Sigma}^{\frac{1}{2}}$的特征值，且特征向量矩阵为$\vect{P},\vect{P}^T\vect{P} = \vect{I}$， 即 $P^T\vect{\Sigma}^{\frac{1}{2}}A\vect{\Sigma}^{\frac{1}{2}}P = diag(\vect{\lambda})$。则定义向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$以及给定一个正数$y$，则二次型随机变量$Y$的概率函数可以通过如下的多项式序列展开进行计算：
	\begin{equation} \label{eq.PROB CDP_infinity}
		p(Y\le y) = \sum_{w=0}^{\infty}(-1)^w c_w \frac{y^{\frac{p}{2} + w}}{\Gamma(\frac{p}{2}+w+1)}, \quad 0<y<\infty,
	\end{equation}
	式中的参数序列$c_w$的计算方法为：
	\begin{equation} \label{eq.PROB_c_w}
		\left\lbrace
		\begin{aligned}
			c_0 &= exp(-\frac{1}{2}\sum_{j=1}^{p}b_j^2)\prod_{j=1}^p(2\lambda_j)^{-\frac{1}{2}} \\
			c_w &= \frac{1}{w}\sum_{r=0}^{w-1}d_{w-r}c_r, \quad w \ge 1
		\end{aligned}
		\right.
	\end{equation}
	以及
	\begin{equation}\label{eq.PROB_d_w}
		d_w = \frac{1}{2} \sum_{j=1}^p (1-wb_j^2)(2\lambda_j)^{-w},\quad w \ge 1
	\end{equation}
	式中，$\Gamma(\cdot)$函数为伽马函数，当输入为正整数时，如$n\in \mathbb{R}^+$，其输出等于$(n-1)!$；$b_j$表示向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$中的第$j$个元素。
\end{theorem}

以及一个关于正整数阶乘计算的定理，
\begin{theorem} \label{th. PROB_factorial_lemma}
	给定一个正整数 $n$, 它的阶乘有如下关系式：
	$$ \left( \frac{n}{\mathit{e}} \right)^n < n! < \mathit{e} \left( \frac{n}{2} \right)^n$$
	式中$\mathit{e}$为自然常数。
\end{theorem} 

\BiSubsection{圆盘观测模型与有限项近似}{The disk-model and finite-term Approximation} 

本小节以给定$k$时刻已知数据预测$k+1$时刻的连通为例，基于通信和观测设备的圆盘模型，将连两个节点之间的连通转变为一个正定二次型随机变量的概率分布，然后使用定理~\ref{th.PROB_PowerSeriesExpansion}计算未来时刻存在连通的概率。由于定理~\ref{th.PROB_PowerSeriesExpansion}中概率的计算公式需要累加无穷多个多项式，因此为了易于计算，仅选用有限个多项式的累加来近似真实概率。本小节最后通过数值仿真的方法分析了有限项近似对概率计算的影响，为后文的对截断误差的理论分析以及基于几何近似的校正方法设计提供了基础。


\begin{figure}[ht]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/diskObservationModel.jpg}
	\bicaption[fig. PROB_DiskModel]{}{圆盘观测模型}{Fig.$\!$}{The Disk Observation Model}
\end{figure}


为了在$k$时刻评估控制输入$\vect{U}^k$，多智能体系统需首先调用合适的估计引擎，得出下一时刻，即$k+1$时刻，系统状态的预测均值向量$\vect{\tilde{X}}^{k+1}$以及协方差矩阵$\vect{\Sigma_{\tilde{X}^{k+1}}^k}$。此时为了预测在$k+1$时刻任意两个节点，例如节点$\vect{p}_1^{k+1}$和节点$\vect{p}_2^{k+1}$，之间通信和观测连通存在的概率，可以从预测值中提取对应节点在$k+1$时刻的预测状态分布。由于本章默认采用高斯噪声假设，因此这两个节点在未来时刻的预测状态也是两个高斯分布函数，即$\vect{p}_1^{k+1}\sim \mathcal{N}(\vect{\tilde{p}}_1^{k+1}, \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}})$以及$\vect{p}_2^{k+1}\sim \mathcal{N}(\vect{\tilde{p}}_2^{k+1}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}})$。因此通过减法定义的上述两个正态分布的差也将符合高斯正态分布函数，即$\Delta \vect{p}_{12}^{k+1} = [\left( \Delta \vect{p}_{12}^{k+1} \right)_x, \left( \Delta \vect{p}_{12}^{k+1} \right)_y]^T \sim \mathcal{N}(\vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}})$。

定义节点$\vect{p}_1^{k+1}$和节点$\vect{p}_2^{k+1}$在未来的$k+1$时刻的真实距离为$$d_{1,2}(k+1) = ||\vect{p}_1^{k+1} - \vect{p}_2^{k+1}||_2$$
则将该距离视作随机变量，它的平方可以表示为距正态分布$\Delta \vect{p}_{12}^{k+1}$的二次型的形式：
$$Y = d_{12}^2(k+1) = \left( \Delta \vect{p}_{12}^{k+1} \right)_x^2 + \left( \Delta \vect{p}_{12}^{k+1} \right)_y^2 = \left( \Delta \vect{p}_{12}^{k+1} \right)^T \vect{A} \Delta \vect{p}_{12}^{k+1}$$
式中矩阵$\vect{A} = \vect{I}$为单位阵。

如图~\ref{fig. PROB_DiskModel}所示，基于通信设备和相对测量设备的圆盘模型，两个节点在未来的$k+1$时刻存在通信和观测的概率，即$s_{12}^{k+1} = 1$的概率，可以等价地转换为
\begin{equation*}
	p\left( s_{1,2}^{k+1} = 1 \right) = p\left( Y\le \rho ^2\right)
\end{equation*}
式中，$\rho$为通信设备和相互测量设备最大的工作半径。

因此对未来时刻连通概率的预测可以转换为两个节点相对距离平方的概率分布，即可以使用定理~\ref{th.PROB_PowerSeriesExpansion}来计算相应的概率。当考虑二维空间的任务规划时，所有节点的状态均为二维向量，因此定理~\ref{th.PROB_PowerSeriesExpansion}中的维度$p=2$。又由于该概率计算公式需要累加无穷多个多项式，在数值计算中是不可行的，因此本章提出仅使用有限项，即前$w = {1,2,...,w_F}$项来做近似计算。由此式~\eqref{eq.PROB CDP_infinity}可变为：
\begin{equation} \label{eq. PROB_CDP_finite}
	p(s_{12}^{k+1} = 1) = p\left( Y \le \rho ^2\right) =  \sum_{w=0}^{w_F}(-1)^w c_w \frac{ （(\rho^2)^{w + 1}}{(w+1)!}
\end{equation}
然而这种通过截取有限项来近似地计算概率的方法必然会引入计算误差。因此上式在应用上的难点是如何保持算法在计算上的可行性的同时减少由有限项近似引入的截断误差。

下面针对有限项近似对概率计算公式\eqref{eq.PROB CDP_infinity}的影响进行讨论，通过仿真实验结果观察到两个主要的问题所在。

\begin{figure}[ht]
	\centering
	\includegraphics[width = \textwidth]{figures_PROB/PerformanceEvaluationTwoNodes1}
	\bicaption[fig. PROB_PerEvaTwoNodes]{}{有限项近似方法计算连通概率存在的问题}{Fig.$\!$}{Indication of the problems for finite-term approximation}\vspace{0em}
\end{figure}

如图~\ref{fig. PROB_PerEvaTwoNodes}所示，本节使用有限项近似的计算结果与基于正态分布抽样点的频率对比，用以侦测有限项近似对频率计算公式~\ref{eq.PROB CDP_infinity}的影响。此处选取两个随机分布节点相对观测的场景。场景一的两个节点正态分布如图~\ref{fig. PROB_PerEvaTwoNodes} A 所示，场景二如图\ref{fig. PROB_PerEvaTwoNodes} D 所示。两个场景中，节点之间的平均距离均为$5m$。两个节点之间的相对距离的分布如图 B 和 E 所示。两个场景的区别在于场景二两个节点之间相对距离的分布比场景一更为稠密。图 C 和 F 为两个节点之间的相对距离的密度分布示意图，两种仿真场景总共抽样了10000个样本，纵坐标为两个节点相对距离与横坐标相等的样本数。 图 G 和图 H 为两个场景中，当观测和通信半径$\rho$从0变化到10的时候，两个节点之间的通信和观测连通存在的真实概率（以10000次抽样的频率为真值参考），以及使用不同$w_F$的有限项近似公式得到的近似概率。

从图~\ref{fig. PROB_PerEvaTwoNodes} G-H中可以观察到两个主要的误差来源。首先是观察子图 E 中当$w_F = 50$和 $w_F = 100$时，有限项近似公式计算得到的近似概率会随着通信和观测设备工作半径$\rho$的增加而产生跳变（$w_F = 50$的变化曲线在约$\rho>6$时突然增加超出概率的$[0,1]$域值，$w_F = 100$的变化曲线在$\rho>8$时突变为负值）。有限项近似的计算方法仅可稳定且准确地计算工作半径$\rho$小于突变点的概率值。突变的方向与使用的多项式个数$w_F$相关。通过观察式~\eqref{eq.PROB CDP_infinity}，该误差的产生源自累加项中存在的$(-1)^k$。该跳变数列的存在使得两个相邻的累加项随设备半径$\rho$的变大而趋向于完全相反增长方向。因此当累加无穷多个多项式计算时，得出的概率可以被限制在合理的$[0 1]$之间。然而当选取有限个多项式的累加作近似计算时，由于忽略了后续的多项式，因此导致朝向一个方向增长的计算数值得不到平衡，因此整体的概率计算向该方向无限增长，从而超出合理范围。

有限项近似导致的第二个计算误差如图~\ref{fig. PROB_PerEvaTwoNodes} H 所示，即当相对距离的分布较为集中，从而导致距离分布的协方差矩阵较小或者近似奇异时，有限项近似方法在临界半径范围内（图 H 中约 $4 < \rho < 6$），无论如何选取累加的多项式数目$w_F$，均无法估计连通存在的概率。
该问题的由来可追溯至定理~\ref{th.PROB_PowerSeriesExpansion}中参数序列$c_w$和向量$\vect{b}$的定义。根据前一节的建模过程以及定理~\ref{th.PROB_PowerSeriesExpansion}，$c_0$的计算需要用$\vect{\Sigma}^{\frac{1}{2}}\vect{A}\vect{\Sigma}^{\frac{1}{2}}$的特征根的倒数。而当二次型随机变量为节点之间相对距离的平方时，$\vect{A} = \vect{I}$，因此$c_0$需要用到$\vect{\Sigma}$特征根的倒数。
另外，向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$也需要用到协方差矩阵$\vect{\Sigma}$的逆。因此当该协方差矩阵的特征值较小或者矩阵近似奇异时，参数序列$c_w$的向量$\vect{b}$的值将被异常放大，从而导致有限项近似的计算偏离正确概率。

\BiSubsection{截断误差稳定性分析}{Stability analysis for truncation error} \label{subsec. ch5sec3.2}

本小节将对有限项近似所引入的截断误差进行理论分析，并给出该误差小于给定阈值的累加项数目选择方法。通过对比定理~\ref{th.PROB_PowerSeriesExpansion}中原始的无穷项累加公式~\ref{eq.PROB CDP_infinity}和有限项近似公式~\ref{eq. PROB_CDP_finite}可以得到截断误差的表达式为：
\begin{equation} \label{eq. PROB_truncationError}
	\begin{aligned}
		\Delta p(s_{1,2}^{k+1}) = \sum_{w=w_F + 1}^{\infty}(-1)^{w} c_w \frac{(y)^{w+1}}{(w+1)!} 
	\end{aligned}
\end{equation}
式中，$y=\rho^2$。

虽然该截断误差也是无穷多个多项式的累加和，然而下面将通过理论推导证明，通过选取合适的$w_F$时，可以任意地减少该截断误差。首先，推论~\ref{cor. PROB_dw} 总结了公式~\eqref{eq.PROB_d_w}定义的参数序列$d_w$的性质如下：
\begin{corollary} \label{cor. PROB_dw}
	给定一个二维正态向量$\vect{X}\sim\mathcal{N}(\vect{\mu},\vect{\Sigma})$，用$\vect{\lambda} = [\lambda_1,\lambda_2]^T$表示协方差矩阵$\vect{\Sigma}$的特征值，矩阵$\vect{P}$为协方差矩阵$\vect{\Sigma}$的特征向量矩阵。定义向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{\frac{1}{2}}\vect{\mu} = [b_1,b_2]^2$，则当下列条件满足时，
	\begin{equation*}
		\left \lbrace
		\begin{aligned}
			&2\lambda_j > 1  \\ 
			&w \ge \frac{1}{b_j^2} + \frac{1}{2\lambda_j - 1} \\
		\end{aligned}	
		\right.
		,\forall j \in \{ 1,2 \}
	\end{equation*}
	公式~\eqref{eq.PROB_d_w}定义的参数序列$d_w$的极限为零，即
	$$\mathop{\lim}_{k\rightarrow \infty} d_k = 0$$
	且序列$d_w$是有界的，即存在正的实数$\bar{d}$,使得
	$$|d_k|<\bar{d},\forall k \in \mathbb{R}^+$$
\end{corollary}

\begin{proof}
	该推论的证明参见附录~\ref{appendix. PROB_dw}。
\end{proof}

然后创建两个新的参数序列$\tilde{c}_w$和$\tilde{d}_w$，其中$\tilde{d}_w$由$d_w$序列的最大值$\bar{d}$组成，即
\begin{equation}\label{eq. PROB_tilde_d_w}
	\tilde{d}_w := [\bar{d}, \bar{d}, ... ,\bar{d}]
\end{equation}

$\tilde{c}_w$的定义分两个部分，首先使$\tilde{c}_w$的初始值与$c_w$的相同，然后使用与公式~\eqref{eq.PROB_c_w}相同的方法计算$\tilde{c}_w$其他的值，即
\begin{equation} \label{eq. PROB_tilde_c_w}
	\begin{aligned}
		\tilde{c}_0 &= c_0 = exp(-\frac{1}{2}\sum_{j=1}^{2}b_j^2)\prod_{j=1}^2(2\lambda_j)^{-\frac{1}{2}} \\
		\tilde{c}_w &= \frac{1}{w}\sum_{r=0}^{w-1} \tilde{d}_{w-r}\tilde{c}_r, \quad w \ge 1
	\end{aligned}
\end{equation}

$\tilde{c}_w$和$\tilde{d}_w$是原始参数序列$c_w,d_w$包络上界。由推论~\ref{cor. PROB_dw}可知，$d_w < \tilde{d}_w, \forall w \in \mathbb{R}^+$，因此$\tilde{d}_w$是$d_w$的包络上界。下面的推论~\ref{cor. PROB_recursive_c_w}给出了$\tilde{c}_w$是$c_w$上界包络的证明，并且推导出了$\tilde{c}_w$的一个递归计算的方法。

\begin{corollary} \label{cor. PROB_recursive_c_w}
	由式~\eqref{eq. PROB_tilde_d_w}和式~\eqref{eq. PROB_tilde_c_w}定义的序列 $\tilde{c}_w$ 的递归计算方程为：
	\begin{equation} \label{eq. PROB_c_w_recursive}
		\tilde{c}_{w+1} = \frac{\bar{d}+w}{w+1} \tilde{c}_{w}
	\end{equation}
	且$\tilde{c}_w$是式~\eqref{eq.PROB_c_w} 定义的参数序列$c_w$的包络上界，即
	$$		|c_k| \le |\tilde{c}_k|, \forall k \in \mathbb{R}^+ $$
\end{corollary}

\begin{proof}
	该推论的证明参见附录~\ref{appendix. PROB_dw}。
\end{proof}

定义$D\in \mathbb{R}^+$为恰好大于$d_w$上界$\bar{d}$的正整数，即$\bar{d} \le D < \bar{d} + 1$。根据$D$的取值，再定义一个新的序列$e_w^D$，

\begin{itemize}
	\item[1):] 当$D=1$时，	
	\begin{equation} \label{eq. PROB_e_w^1}
		e^D_w = 	\frac{c_w}{(w+1) }, \forall w \in \mathbb{R}^+
	\end{equation}
	
	\item[2):] 当$D \ge 2$时，
	\begin{equation} \label{eq. new series e_k^D, D>2}
		e^D_w = \left\lbrace 
		\begin{array}{ll}
			c_w, & \forall w \le (D - 2) \\
			\frac{c_w}{\prod_{j=1}^{D} (w+2-j) }, & \forall w \ge (D - 1)
		\end{array}	\right.
	\end{equation}
\end{itemize}

则一个关于序列$e_w^D$的性质总结为如下推论：
\begin{corollary} \label{cor. PROB_e_w^D}
	假设已知正整数$D$，则序列$e_w^D$的极限是有界的，即
	$$\mathop{\lim}_{w \rightarrow \infty} |e_w^D| = \mathit{Const.} ,$$
	且$e_w^D$的最大值满足$|e_w^D| \le \tilde{c}_{D-1}, \forall w \in \mathbb{R}^+$
\end{corollary}

\begin{proof}
	该推论的证明参见附录~\ref{appendix. PROB_e_w^D}。
\end{proof}

由以上三个推论的支撑，式~\eqref{eq. PROB_truncationError}所定义的截断误差满足如下定理:

\begin{theorem} \label{th. PROB_adaptive_main_theorem}
	给定一个二维正态向量$\vect{X} \sim \mathcal{N}(\vect{\mu},\vect{\Sigma})$，用$\vect{\lambda} = [\lambda_1,\lambda_2]^T$表示协方差矩阵$\vect{\Sigma}$的特征值，矩阵$\vect{P}$为协方差矩阵$\vect{\Sigma}$的特征向量矩阵。定义向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{\frac{1}{2}}\vect{\mu} = [b_1,b_2]^T$，并且按照式~\eqref{eq.PROB_c_w}定义参数序列$c_w$，以及按照式~\eqref{eq.PROB_d_w}定义$d_w$，令$D \in \mathbb{R}^+$表示$d_w$上界的正整数。首先假设
	$$2\lambda_j >1,\forall j \in \{1,2\} $$
	则当任意给定一个误差阈值$\delta > 0$，由式~\eqref{eq. PROB_truncationError}所定义的截断误差$\Delta p$ 可以被减少到该误差阈值之下，即
	\begin{equation*}
		\Delta p(s_{1,2}^{k+1} = 1) < \delta
	\end{equation*}
	仅当有限项近似中的累加项数目$w_F$满足如下条件时：
	\begin{equation}  \label{eq. PROB_w_F_conditions}
		\left\lbrace
		\begin{aligned}
			&w_F \ge \frac{1}{b_j^2} + \frac{1}{2\lambda_j - 1} \forall j \in \{ 1,2 \} \\
			&w_F \ge y(\mathit{e})^2 + D - 2 \\
			& w_F \ge 3D - ln \left( \frac{\delta}{c_0y^D} \right) - 2	 	
		\end{aligned}
		\right.		
	\end{equation}
	式中$\mathit{e}$表示自然常数。
\end{theorem}

\begin{proof}
	基于推论~\ref{cor. PROB_dw},~\ref{cor. PROB_recursive_c_w},~\ref{cor. PROB_e_w^D}，则式~\eqref{eq. PROB_truncationError}定义的截断误差可作如下展开：
	\begin{equation*}
		\begin{aligned}
			\Delta &p(s_{1,2}^{k+1}) = \sum_{w=w_F+1}^{\infty}(-1)^w c_w \frac{y^{w+1}}{(w+1)!} \\
			&\le  \sum_{w=w_F+1}^{\infty} \frac{|c_w|}{\prod_{j=1}^{D}(w+2-j)} \frac{y^{w+1}}{(w+1-D)!} \\
			&\le \sum_{w=w_F+1}^{\infty} |\tilde{e}^D_w|\frac{y^{w+1}}{(w+1-D)!} \\
			&\le \tilde{c}_{D-1} \frac{y^{w_F+2}}{(w_F+2-D)!} \left(1 + \sum_{w=1}^{\infty} \frac{y^{w}}{\prod_{j=1}^{w} (w_F + 2+ j -D) }  \right) \\
			&< \tilde{c}_{D-1} \frac{y^{w_F+2}}{(w_F+2-D)!} \left(1 + \sum_{w=1}^{\infty} \frac{y^{w}}{ (w_F -D)^w }  \right)
		\end{aligned}
	\end{equation*}
	最后一行不等式右侧的无穷累加部分可以通过巴塞尔问题来求解，即寻找合适的$w_F$使
	$$		\frac{1}{ \left( \frac{w_F-D}{y} \right)^k } \le \frac{1}{k^2}$$
	则根据巴塞尔问题可得，
	$$		\left(1 + \sum_{w=1}^{\infty} \frac{y^{w}}{ (w_F -D)^w }  \right) \le \left( 1 + \sum_{w=1}^{\infty} \frac{1}{w^2} \right)  = 1 + \frac{\pi^2}{6} $$
	因此使得该无穷累计项满足巴塞尔问题的条件为，
	\begin{equation*}
		\begin{aligned}
			&\frac{1}{ \left( \frac{w_F-D}{y} \right)^w } \le \frac{1}{w^2}\\
			\Rightarrow\quad & \left( \frac{w_F-D}{y} \right)^w \ge w^2 \\
			\Rightarrow\quad & \ln(w_F-D) - \ln(y) \ge \frac{2\ln(w)}{w}. 	
		\end{aligned}
	\end{equation*}
	由于函数$\frac{\ln(w)}{w}$的最大值为$\frac{1}{\mathit{e}}$， 由此可得
	\begin{equation} \label{eq. PROB_w_F_Basel_condition}
		\begin{aligned}
			&\ln(w_F - D) - \ln(y) \ge \frac{2}{e} \\
			\Rightarrow \quad & w_F \ge e^{\frac{2}{e}}y + D
		\end{aligned}
	\end{equation}
	
	\begin{table}[tbp]
		\bicaption[tab. PROB_g_w_F]{}{不同条件下$g(w_F)$低于四组域值时所需的$w_F$值}{Table$\!$}{Values of $w_F$ when $g(w_F)$ is less than a certain threshold given finite $y$ and $D$}
		\vspace{0.5em}
		\centering
		\wuhao		
		\begin{tabular}{ccccccc}
			\toprule[1.5pt]
			\multirow{2}{*}{条件} & \multicolumn{4}{c}{ $g(w_F) \le$ } \\ \cline{2-5}
			& $10^{-3}$ & $10^{-5}$ & $10^{-10}$ & $10^{-20}$ \\ 
			\midrule[1pt]
			y = 25($\rho$ = 5), D = 10               & 108   & 109   & 118    & 133    \\
			y = 25($\rho$ = 5), D = 20              & 139   & 141   & 149    & 163    \\
			y = 100($\rho$ = 10), D = 10              & 327  & 327  & 340   & 359   \\
			y = 100($\rho$ = 10), D = 20             & 375  & 378  & 387   & 405   \\ 
			\bottomrule[1.5pt]
		\end{tabular}
	\end{table}
	
	下面考虑截断误差展开的不等式的中间项，令$g(w_F) = \frac{y^{w_F+2}}{(w_F+2-D)!}w_F$，则将该序列相邻两项元素做差可得，
	$$
	\begin{aligned}
		& g(w_F+1) - g(w_F)\\ 
		&= \frac{y^{w_F+3}}{(k_m+3-D)!} - \frac{y^{w_F+2}}{(w_F+2-D)!}\\
		& = \frac{y^{w_F+2}}{(w_F+2-D)!} \left( \frac{y}{w_F+3-D} - 1\right)
	\end{aligned}
	$$
	因此，当$y=\rho^2$和$D$全都已知时，令$w_F>y+D-3$即可让$g(w_F+1)<g(w_F)$，因此$w_F$的取值越大，$g(w_F)$就越小。$g(w_F)$随$w_F$的增加而变小的趋势如表~\ref{tab. PROB_g_w_F}所示。
	
	在满足巴塞尔条件~\eqref{eq. PROB_w_F_Basel_condition}后，截断误差为
	\begin{equation} \label{eq. PROB_trunctionErrorAfterBasel}
		\begin{aligned}
			\Delta p(s_{1,2}^{k+1}) 
			&< \tilde{c}_{D-1} \left( 1+\frac{\pi^2}{6} \right)g(w_F)
		\end{aligned}
	\end{equation}
	由推论~\ref{cor. PROB_recursive_c_w} 可知，$\tilde{c}_{D-1}$是有限值，因此经过上述的定性分析之后可得，增加$w_F$的取值可以减少截断误差。
	
	下面将定量地分析$w_F$的取值如何影响截断误差，即确定合适的$w_F$，当给定一个截断误差阈值的$\delta$时，使得$\Delta p(s_{1,2}^{k+1}) < \delta$。
	
	首先分析$\tilde{c}_{D-1}$，根据推论~\ref{cor. PROB_recursive_c_w}可知
	\begin{equation*}
		\begin{aligned}
			\tilde{c}_{w} & = \frac{w-1+\bar{d}}{w}\tilde{c}_{w-1} \\
			& = \frac{(w-1+\bar{d})(w-2+\bar{d})}{w(w-1)}\tilde{c}_{w-1} \\
			&= \frac{\prod_{i=0}^{w-1} (\bar{d}+i) }{\prod_{j=0}^{w-1} (j+1) } \tilde{c}_0	
		\end{aligned}
	\end{equation*}
	由于$\bar{d} < D$，$ \prod_{i=0}^{w-1} (\bar{d}+i) \le \prod_{i=0}^{w-1} (D+i) = \frac{(D+w-1)!}{(D-1)!}$，因此可知
	\begin{equation*}
		\tilde{c}_{w} \le \frac{(D+w-1)!}{w!(D-1)!} \tilde{c}_0, \forall w \in \mathbb{R}^+
	\end{equation*}
	因此$\tilde{c}_{D-1}$的上界为
	\begin{equation*}
		\tilde{c}_{D-1} \le \frac{(2D-2)!}{(D-1)!(D-1)!} \tilde{c}_0
	\end{equation*}	
	由定理~\ref{th. PROB_factorial_lemma}可得
	\begin{equation} \label{eq. PROB_maximum_tile_c_D_1}
		\begin{aligned}
			\tilde{c}_{D-1} &\le \frac{(2D-2)!}{(D-1)!(D-1)!} \tilde{c}_0 \\
			& \le \frac{e(D-1)^{2D-2}}{ \left( \frac{D-1}{e} \right)^{D-1} \left( \frac{D-1}{e} \right)^{D-1} } \tilde{c}_0 \\
			& \le e^{2D-1} \tilde{c}_0 = e^{2D-1} c_0
		\end{aligned}
	\end{equation}	
	
	下面分析$g(w_F)$与$w_F$取值的具体关系。使用变量替换$\alpha = w_F + 2 - D$，则再次引用定理~\ref{th. PROB_factorial_lemma}可得：
	\begin{equation} \label{eq. PROB_g(w_F)_alpha}
		\begin{aligned}
			& g(w_F) = \frac{y^{\alpha+D}}{\alpha!} \\
			& \le \frac{y^{\alpha+D}e^{\alpha}}{{\alpha}^{\alpha}} = \left( \frac{ye}{\alpha} \right)^{\alpha} y^D
		\end{aligned}		
	\end{equation}
	将公式~\eqref{eq. PROB_g(w_F)_alpha}和公式~\eqref{eq. PROB_maximum_tile_c_D_1}代入公式~\eqref{eq. PROB_trunctionErrorAfterBasel}，则当给定截断误差的阈值$\delta$时，为了使$\Delta p(s_{1,2}^{k+1}) < \delta$，需保证
	\begin{equation*}
		\Delta p(s_{1,2}^{k+1}) \le e^{2D-1}c_0 (1+\frac{\pi^2}{6}) \left( \frac{ye}{\alpha} \right)^{\alpha} y^D < \delta
	\end{equation*} 
	由于$1+\frac{\pi^2}{6} < \mathit{e}$，因此上式可转换为：
	\begin{equation*}
		\begin{aligned}
			e^{2D} \left( \frac{ye}{\alpha} \right)^{\alpha} < \frac{\delta f}{c_0 y^D}
		\end{aligned}
	\end{equation*}
	
	首先令$\alpha \ge e^2y \Rightarrow w_F \ge e^2y + D - 2$,此时
	\begin{equation*}
		e^{2D} \left(\frac{ye}{\alpha} \right)^\alpha \le e^{(2D-\alpha)}
	\end{equation*}
	然后令
	\begin{equation*}
		\begin{aligned}
			& e^{(2D-\alpha)} < \frac{\delta f}{c_0y^D} \\
			\Rightarrow & \alpha > 2D - ln \left( \frac{\delta f}{c_0y^D} \right) \\
			\Rightarrow & w_F > 3D - ln \left( \frac{\delta f}{c_0y^D} \right) - 2
		\end{aligned}
	\end{equation*}
	
	因此，综合巴塞尔问题对$w_F$的要求，可以总结出令$\Delta p(s_{1,2}^{k+1}) < \delta$时$w_F$的选取规则，如式~\eqref{eq. PROB_w_F_conditions}所示。至此，定理~\ref{th. PROB_adaptive_main_theorem}得证。
\end{proof}

\BiSubsection{基于平移膨胀的协方差校正}{Translationally approximated covariance expansion }

定理~\ref{th. PROB_adaptive_main_theorem}表明，基于有限项近似的概率计算方法在选取合适的累加项个数$w_F$后，其计算结果可以任意地接近真值。然而以上的分析全部基于协方差矩阵特征值大于$\frac{1}{2}$的假设。在实际的多智能体协同任务中，智能体对未来状态的预测是实时改变的，因此协方差矩阵特征值也是实时变化的。在多智能体系统的协同任务中，处于临界状态的连通往往决定了系统的最佳性能表现。因此关注临界的通信和工作半径状态下节点之间连通的计算过程是十分重要的。由以上分析可知，有限项近似的概率计算方法在相对距离分布的协方差矩阵较小或者近似奇异时表现不佳。因此，解决该问题的思路是将较小或者近似奇异的矩阵进行等比例膨胀，将膨胀后正常的矩阵带入有限项近似的概率计算公式，作近似计算。

本节提出一种基于平移膨胀的协方差校正方法，下面将详细介绍协方差近似膨胀方法的设计过程。如图~\ref{fig. PROB_Trace}所示，红色实线的椭圆$E_{r1}$代表两个节点初始情况下相对距离分布的协方差矩阵$\vect{\Sigma}$的$3\delta$置信域。此时两个节点上配备的通信和测量设备的工作半径为$\rho_1$，则以坐标原点为中心，$\rho_1$为半径的通信测量圆与$E_{r1}$相交于$A$和$C$点。假设$E_{r1}$内部存在一个小的椭圆$E_{b1}$与半径为$\rho_1$通信测量圆相切于$B$点。由于$E_{r1}$所代表的协方差矩阵较小或者是近似奇异，因此将$E_{r1},E_{b1}$沿着$OO_1$共线方向等比例放大$\beta$倍，获得新的置信域$E_{r2}$，内部椭圆$E_{b2}$以及对应的交点$F,D$，相切点$E$。通过等比例放大，新的置信椭圆$E_{r2}$对应的协方差矩阵不会存在奇异现象。

则由相关定义可知，两个节点之间存在通信和观测连通的概率等于节点之间相对距离小于设备工作半径，即图中通信观测圆与置信椭圆之间的交集部分，其数学表述为：
\begin{equation*}
	p(s_{12}^{k+1} = 1) \approx \frac{S(\widearc{ABC} \cap E_{r1})}{S(E_{r1})} 
\end{equation*} 
式中$S(\cdot)$为计算输入参数的面积。
根据相关的几何知识易得，图中$\{O,O_1,O_2\},\{O,A,D\},\{O,C,F\},\{O,B,E\}$等四组点群，每组点群内的点都处在同一直线上。则根据几何相似可知，
\begin{equation}
	\frac{S(\widearc{DEF} \cap E_{r2})}{S(E_{r2})} =  \frac{S(\widearc{ABC} \cap E_{r1})}{S(E_{r1})}
\end{equation} 
因此，基于平移膨胀的协方差校正方法可以使用膨胀后新的且满足定理~\ref{th. PROB_adaptive_main_theorem}的增广协方差矩阵进行有限项近似的概率计算，此时的等效半径$\rho_{2}$将大于最初的工作半径$\rho_{1}$。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/Trace}
	\bicaption[fig. PROB_Trace]{}{平移膨胀的协方差校正方法}{Fig.$\!$}{Indication of the translationally approximated covariance expansion}\vspace{0em}
\end{figure}


\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwResult{ $\vect{\mu}_2,\vect{\Sigma}_2,\rho_2$ \; }
	\KwData{ $\vect{\mu}_1,\vect{\Sigma}_1,\rho_1$ \;}
	
	$\lambda_{min} \leftarrow$ 计算 $\vect{\Sigma}_1$ 的最小特征根  \;
	
	\eIf{$\lambda_{min} < \frac{1}{2}$} {
		$\beta = \frac{1}{\lambda_{min}},$
		$\vect{\Sigma}_2 = \beta \vect{\Sigma}_1,$
		$\vect{\mu}_2 = \sqrt{\beta} \vect{\mu}_1,$ 
		$\rho_2 = \sqrt{\beta}  \rho_{1}$ \;
	}
	{		$\vect{\Sigma}_2 =  \vect{\Sigma}_1,$
		$\vect{\mu}_2 =  \vect{\mu}_1$ 
		$\rho_2 =  \rho_{1}$ \;
	}
	
	
	\textbf{Return} $\vect{\mu}_2,\vect{\Sigma}_2,\rho_2$
	\AlgoBiCaption{平移协方差膨胀校正算法}{The algorithm of translationally approximated covariance expansion}
	\label{alg. Trace}
\end{algorithm}

本节所提出的协方差平移膨胀是基于高斯噪声假设。只有在高斯噪声条件下才可以使用椭圆置信区间来表示协方差矩阵。因为在高斯模型下理论上将会有超过$99.7\%$的抽样点落在置信椭圆内。

基于平移膨胀的方法也可以被反向使用，即当通信和观测半径$\rho_{1}$比较大，但是协方差矩阵的分布较为分散时，可以反向地使用该膨胀方法，把一个分散的协方差矩阵缩小但是依然符合定理~\ref{th. PROB_adaptive_main_theorem}对协方差矩阵的要求，但是可以换取一个较小的等效半径$\rho_{2}$。该方法的好处是由于$y=\rho^2$在有限项近似的计算公式中被频繁使用，一个较大的工作半径$\rho_{1}$会积累数值计算误差，从而导致最终的结果偏离正确概率。因此较小的$\rho_{2}$将减轻此类问题。

下面采用与协方差近似膨胀相同的方法计算膨胀系数$\beta$，即确保协方差矩阵$\vect{\Sigma}$的最小特征值不小于1。下面将该算法的伪码总结为算法~\ref{alg. Trace}。

\begin{algorithm*}[htbp]
	\SetAlgoLined
	\KwIn{ 	$\qquad p(s_{1,2}^{k+1}=1)$  \;}
	\KwOut{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}, \rho_{1}, \delta$ \;}
	
	\SetKwProg{Fn}{Function}{:}{end}	
	\Fn{ APSE( \textbf{Input} ) } {
	
		$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
		$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
		
		$[\rho_{3\delta}^{min} \rho_{3\delta}^{max}] \leftarrow 3\delta$ 置信域计算( $\vect{\Sigma}_{\Delta}$ )
		
		\uIf{$\rho < \rho_{3\delta}^{min}$}{
			\textbf{Return} $p(Y<\rho^2) = 0$ \;
		} 
		\ElseIf{$\rho > \rho_{3\delta}^{max}$}{
			\textbf{Return} $p(Y<\rho^2) = 1$ \;
		} 
		$\left[ \vect{\mu},\vect{\Sigma},\rho_2 \right] \leftarrow \text{调用算法~\ref{alg. Trace}对协方差矩阵进行校正}(\vect{\mu}_{\Delta},\vect{\Sigma}_{\Delta},\rho_1)$ \;
		
		$\vect{\lambda} = [ \lambda_1, \lambda_2 ]^T \leftarrow$ 计算协方差矩阵 $\vect{\Sigma}$的特征值 \;
		
		$\vect{P} \leftarrow$ 计算协方差矩阵$\vect{\Sigma}$的特征向量矩阵 \;
		$\vect{b} = [b_1,b_2]^T = P^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$, $y = \rho^2 $, $w_F = 50$ \tcc*{赋予$w_F$初始值} 
		
		\For{$r=0;r\le w_F;r=r+1$}{
			$d_r \leftarrow $ 调用公式~\eqref{eq.PROB_d_w}, $\forall r = 1,...,w_F$\;
			$c_r \leftarrow $ 使用 $d_r$ 并调用公式~\eqref{eq.PROB_c_w} $\forall r = 0,...,w_F$\;
			$y_r = (-1)^r \dfrac{y^{r+1}}{\Gamma(r+2)}, r = 0,...,w_F $  \;
			\If{$d_r$达到最大值}{
				$D = d_r$ \tcc*{该条件语句仅进入并执行一次}
				$w_{F1} = \frac{1}{b_1^2} + \frac{1}{2\lambda_1 - 1}$,$w_{F2} = \frac{1}{b_2^2} + \frac{1}{2\lambda_2 - 1}$ \;
				$w_{F3} = e^2y + D-1$, $w_{F4} = 3D - ln \left( \frac{\delta}{c_0y^D} \right) - 2$ \;
				
				$w_F = \text{max}(k_{m1},k_{m2},k_{m3},k_{m4},k_m)$ \;
			}
		}
		$ p(s_{1,2}^{k+1}=1) \leftarrow $ 基于 $c_r, y_r$调用公式~\eqref{eq. PROB_CDP_finite}预测连通概率 \;
		
		\textbf{Return} $p(s_{1,2}^{k+1}=1)$
	}
	\AlgoBiCaption{APSE：基于自适应参数选取的有限项近似概率改进算法}{APSE: Adaptive Power Seris Extension algorithm for the prediction of probability}
	\label{alg. PROB_APSE}
\end{algorithm*}

\BiSubsection{基于自适应参数选择的连通概率预测算法}{Probability prediction algorithm: the adaptive power series extension}

前两节给出了两种有限项近似的改进策略。首先是通过严格的理论推导，得到了有限项近似方法累加项数目$w_F$的自适应选取策略。然后是针对自适应策略中存在的对协方差矩阵的约束，设计了一种改进的协方差膨胀方法。基于以上分析，本节将基于自适应参数选取的改进概率预测算法的伪码总结为算法~\ref{alg. PROB_APSE}。

按照从上到下的顺序，算法~\ref{alg. PROB_APSE}首先需要两个节点对未来$k+1$时刻的预测状态$\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$以及节点上装配的相对通信和测量设备的工作距离$\rho$和期望的截断误差精度$\delta$作为算法输入。最终返回两个节点在未来时刻存在通信和相对测量连通的概率，即$p(s_{1,2}^{k+1})$。算法的$3-9$行是一种基于高斯概率性质的$3\delta$区间校正方法。基于高斯噪声假设，两个节点的分布$\vect{p}_1^{k+1}$和$\vect{p}_2^{k+1}$均为正态分布，因此两个节点之间通过减法定义的矢量$\Delta \vect{p}_{12}^{k+1}$也是一个二维正态分布的向量。如图~\ref{fig. PROB_3delta}所示，根据高斯正态分布函数的性质可知，$\Delta \vect{p}_{12}^{k+1}$分布中$99.7\%$的抽样点均将落入其$3\delta$置信区间中。因此，此处的$3 \delta$区间校正方法是指在获取节点之间的减法矢量$\Delta \vect{p}_{12}^{k+1}$的分布后，通过提前计算该分布的$3\delta$区间，得到边界$[\rho_1,\rho_2]$，则当节点的通信观测设备的工作半径$\rho<\rho_1$时，将节点之间存在连通的概率直接赋值为0，否则当$\rho>\rho_2$时，将概率直接赋值为1。使用该方法的目的是为了在节点之间的相对距离远大于或远小于设备的工作半径时，直接判定连通度，从而减少算法的计算复杂度。
\begin{figure}[ht]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/3_delta_figure}
	\bicaption[fig. PROB_3delta]{}{$3\delta$修正方法示意图}{Fig.$\!$}{Indication for the $3\delta$ regulation }\vspace{0em}
\end{figure}

然后第$10$行调用了改进的平移膨胀协方差校正算法，以确保节点之间相对距离分布的协方差矩阵符合定理~\ref{th. PROB_adaptive_main_theorem}的要求。$11-13$行是有限项近似公式~\eqref{eq. PROB_CDP_finite}的前期变量处理过程。$14-26$行执行有限项近似的概率计算过程，其中$19-24$行为有限项近似中累加项数目$w_F$的自适应选择方法，对应于定理~\ref{th. PROB_adaptive_main_theorem}中的选取条件。在自适应地选取$w_F$的值之前，算法为$w_F$选定了一个初始值，即$14$行中的$w_F=50$。该初始化过程实际上不是必需的，仅是为了初始情况下方便程序对内存进行分配。另外，该算法假设在给定的初始累加项数目里（此算法是50项），$d_w$能够寻找到该序列的最大值，因此可以进入$19-24$行的条件语句。由于最终的$w_F$的取值是将定理~\ref{th. PROB_adaptive_main_theorem}中所有的候选值与初始值作比较并取最大值，因此有限项近似的计算公式中将至少累加50项

\BiSubsection{概率预测算法精准度验证仿真}{Simulation Results for the Accuracy of Predicted Probability} \label{subsec. ch5sec5.3.5}
本节测试算法~\ref{alg. PROB_APSE}对节点之间连通存在概率预测的性能表现，此处考虑两类性能指标，即算法对连通概率预测的准确性和算法的数值运算复杂度。为了验证算法的有效性，此处的仿真验证将算法~\ref{alg. PROB_APSE}与采取不同抽样次数的随机抽样方法作对比。随机抽样方法是一种通过拟合样本分布，经过重复抽样，统计事件发生的频率来模拟概率的方法。将使用随机抽样方法预测连通概率的过程总结为算法~\ref{alg. RandomSampleDistributionAlg}。

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwOut{ $p(s_{1,2}^{k+1}=1)$  }
	\KwIn{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}},\rho, dim$ \;}
	
	
	\SetKwProg{Fn}{Function}{:}{end}	
	\Fn{ RandomSampling( \textbf{Input} ) } {
		$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
		$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
		
		Count = 0 \;
		\While{$k \le dim$}{
			$\vect{\mu}$ $\leftarrow$ 根据$\mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$的分布进行随机抽样 \;
			\If{ $||\vect{\mu}||_2$ $ \le \rho$}
			{Count = Count + 1 \;}
		}
		\textbf{Return} $p(s_{1,2}^{k+1}=1) = \frac{Count}{dim}$\;
	}
	\AlgoBiCaption{随机抽样算法}{Randomly sampling algorithm}
	\label{alg. RandomSampleDistributionAlg}
\end{algorithm}

算法~\ref{alg. PROB_APSE}和算法~\ref{alg. RandomSampleDistributionAlg}的仿真实现是基于 Windows 系统下的 MATLAB R2019a 软件，运行于i7-9750H@2.6GHz 处理器和8GB内存的计算平台。由于算法~\ref{alg. PROB_APSE}在具体的程序实现上是运动MATLAB的\textit{for-loop}语句，而算法~\ref{alg. RandomSampleDistributionAlg}的实现是采用MATLAB软件的内置函数\textit{mvrnd}函数，因此算法~\ref{alg. PROB_APSE}在数值计算效率上存在算法层面和优化层面的提升空间。

下面将仿真的具体设置介绍如下。测试两个节点$A$和$B$，它们的真实位置固定为$\vect{p_A} = [0.5,1]^T$ 以及 $\vect{p}_B=[2,2.5]^T$，因此节点$A，B$之间的真实距离为$||\vect{p}_A - \vect{p}_B||_2 \approx 2.21$。然后随机生成两个正定对称矩阵$\vect{P}_A>0,\vect{P}_B>0$作为对$A,B$两个节点估计位置分布的协方差矩阵。对任意的一组协方差分布，设备的工作半径共测试$\rho$为$\rho \in \vect{S}_{\rho}:=\{ 0.1:0.1:6 \}$集合中的60个数值，以此验证算法对连通存在的确定存在状态、确定不存在状态以及临界状态的表现性能。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.9\textwidth]{figures_PROB/32RealDistributionAndErrorsSingleCase}
	\bicaption[fig. PROB_RealDistributionAndErrorsSingleCase]{}{典型案列的真实分布}{Fig.$\!$}{Real distribution of the representative case}\vspace{0em}
\end{figure}

\begin{figure}[htbp]
	\centering	
	\subfigure{\label{subfig. discretizationArea}}\addtocounter{subfigure}{-2}
	\subfigure[RMSE for predicted probability error]{\subfigure[概率预测误差的RMSE]{\includegraphics[width = 0.45\textwidth]{figures_PROB/32RMSESingleCase1}}}	
	\null\hfill
	\subfigure{\label{subfig. discretizationBoundary}}\addtocounter{subfigure}{-2}
	\subfigure[Run time]{\subfigure[数值计算复杂度]{\includegraphics[width = 0.45\textwidth]{figures_PROB/32RunningTimeSingleCase1}}}
	
	\bicaption[fig. PROB_performanceRepresentativeCase]{}{典型案例下算法的性能表现与对比}{Fig.$\!$}{Performance comparison for the representative case}
	\vspace{-0em}
\end{figure}

然后对于任意的$\rho \in \vect{S}_{\rho}$，将$\vect{p}_A,\vect{p}_B,\vect{P}_A,\vect{P}_B$和$\rho$作为算法~\ref{alg. PROB_APSE}以及算法~\ref{alg. RandomSampleDistributionAlg}的输入。对于算法~\ref{alg. PROB_APSE}最后的待输入参数，截断误差的期望阈值$\delta$，此处设置为固定值$10^{-10}$。作为对照组，随机抽样算法~\ref{alg. RandomSampleDistributionAlg}共测试六组不同的抽样值，即$dim=[10,10^2,10^3,10^4,10^5,10^6]$。概率的真实值使用$dim=10^7$的算法~\ref{alg. RandomSampleDistributionAlg}得出的频率值来近似。


首先测试一组典型案列，随机生成两个协方差矩阵为
$$\vect{P}_A = \left[ \begin{matrix}
	0.371629432909261&-0.100755606016948 \\
	-0.100755606016948&	0.263717689138266
\end{matrix} \right],$$
以及
$$\vect{P}_B = \left[ \begin{matrix}
	0.196735055559263&0.370216978853931 \\
	0.370216978853931&0.720201443148860
\end{matrix} \right].$$
因此节点$A$和$B$之间的相对距离分布为$\Delta \vect{p} \sim \mathcal{N}(\vect{p_A} - \vect{p}_B, \vect{P}_A+\vect{P}_B)$。通过计算可知，协方差$\vect{P}_A+\vect{P}_B$的最小特征值为 $0.435875962912945<0.5$，因此，该典型案例在调用算法~\ref{alg. PROB_APSE}时将激活算法~\ref{alg. Trace}进行协方差膨胀的校正过程。

典型案例中两个节点的分布如图~\ref{fig. PROB_RealDistributionAndErrorsSingleCase} A 所示，两个节点位置的相对分布如图中 B 所示。子图 C 展示了针对典型案例，实验组和6个对照组的概率预测方法所得到的概率预测结果，与真实概率之间之间的误差，随设备工作半径$\rho$变化的分布。图~\ref{fig. PROB_performanceRepresentativeCase} 是对实验组和6个对照组在综合了集合$\vect{S}_{\rho}$中60个半径$\rho$取值后，获得的概率预测误差的RMSE值以及每组方法的平均运算时间（包括真实概率的仿真组）。

\begin{figure}[htbp]
	\centering	
	\subfigure{\label{subfig. discretizationArea}}\addtocounter{subfigure}{-2}
	\subfigure[Distribution of RMSE for predicted probability error]{\subfigure[概率误差RMSE分布]{\includegraphics[width = 0.45\textwidth]{figures_PROB/81RMSEALLCases1}}}	
	\null\hfill
	\subfigure{\label{subfig. discretizationBoundary}}\addtocounter{subfigure}{-2}
	\subfigure[Distribution of the run time]{\subfigure[数值计算复杂度分布]{\includegraphics[width = 0.45\textwidth]{figures_PROB/81RunningTimeALLCases1}}}
	
	\bicaption[fig. PROB_performance200Distribution]{}{200次打靶仿真中算法的概率预测表现性能分布}{Fig.$\!$}{Performance distribution of prediction algorithm over 200 repeated trials}
\end{figure}

然后重复上述典型案例的仿真过程，进行200打靶，每次随机生成节点的位置状态协方差矩阵$\vect{P}_A,\vect{P}_B$。记录两种性能表现，分别为平均概率预测误差的RMSE值以及包括真实概率组在内的每种方法的平均运算时间。此两种性能表现在200次的重复实验中的分布如图~\ref{fig. PROB_performance200Distribution}。由仿真结果可得，本节设计的基于自适应参数选取的有限项近似的概率预测算法~\ref{alg. PROB_APSE}在概率预测精度上优于所有的对照组。其数值计算复杂度上介于使用$10^{3}$-$10^4$个抽样的随机抽样算法的计算复杂度之间。

\BiSection{基于拓扑概率预测的构型规划算法}{Configuration planning for CL using an exact network probabiltiy} \label{subsec. ch5sec4}
前一节详细介绍了一种圆盘观测模型连通概率预测方法的设计过程，总结了算法~\ref{alg. PROB_APSE}，并通过数值仿真验证了该算法在数值计算复杂度和连通概率预测精度上的优越性。以前一章所设计的构型规划框架为基础，本节将研究把拓扑概率预测算法应用于构型规划问题，通过精确地预测未来规划期限内多智能体系统内部观测拓扑，提升规划框架中目标函数计算的准确性，从而进一步提升构型规划方法的性能。

依据目标函数在式~\eqref{eq. ACTCO_objective function2}中的分解，当给定待评估的控制输入以及网络的未来观测拓扑时，目标函数的具体计算过程可总结为算法~\ref{alg. EvaluateObjectiveFucntionProbCon}所示。

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwIn{$\rho$传感器最大测量范围； $\delta$ 连通概率计算误差阈值 \\
		\quad $b(\vect{X}^{k}) \sim \mathcal{N}(\vect{\bar{X}}^{k}, \vect{\Sigma}^{k})$:  多智能体系统在规划时刻$k$的信念分布$b(\vect{X}^{k})$\\ 
		\quad $b(\vect{X}^{0}) \sim \mathcal{N}(\vect{\bar{X}}^{0}, \vect{\Sigma}^{0})$: 多智能体的初始信念 $b(\vect{X}^{0})$  \\
		\quad $\mathcal{H}^{0:k+L|k} = \{ \mathcal{H}^{0:k}, \vect{U}^{k:k+L-1} \}$ : 包含待评估控制输入的历史数据集合  \\
		\quad $\mathcal{E}^{k+1:k+L}_{iter}$: 给定的网络拓扑;\\
		%	\quad $r_L$: maximum loop times for the GN iteration }
	}
	\KwOut{\quad $J_k^{iter}$: 给定$\mathcal{E}^{k+1:k+L}_{iter}$和$\vect{U}^{k:k+L-1}$的目标函数值 }
	
	\SetKwProg{Fn}{Function}{:}{end}	
	\Fn{ EvaluateObjectiveProbCon( \textbf{Input} ) } {
		依据最大似然假设预测实际观测数据 $\vect{\vec{Z}}^{k+1:k+L}_{ML} \leftarrow $ 式~\eqref{eq. ACOTO_MLobservationPrediction1}  \;
		从输入数据集$\mathcal{H}^{0:k}$中提取历史观测数据 $\vect{Z}^{0:k}$ \;
		$\vect{Z}^{0:k+L}_{ML} =\{ \vect{\vec{Z}}^{k+1:k+L}_{ML}, \mathcal{E}^{k+1:k+L}_{iter}  \}  \bigcup \vect{Z}^{0:k}$ \;
		\uIf(\tcc*[f]{重复调用 \text{\textbf{算法~\ref{alg. FilteringEngineEstimation}}} }){Filtering Engine}{
			\For{$l=0:L-1$}{
				$(\vect{\bar{X}}^{k+l+1}, \vect{\Sigma}^{k+l+1}) = FilteringEngineEKF(\vect{\bar{X}}^{k+l},\vect{\Sigma}^{k+l},\vect{U}^{k+l},\vect{Z}^{k+l+1}_{ML}); $ 	\\
				计算先验估计$b(\vect{X}^{k+l+1|k+l}) \leftarrow$ 式~\eqref{eq. prioriUpdateEKF}\;
			}
			$b(\vect{{X}}^{k+1:k+L|k}) = \{ b(\vect{X}^{k+i+1|k+i}) | \forall 0 \le i \le L-1 \}$ \;
		}
		\ElseIf(\tcc*[f]{调用 \text{\textbf{算法~\ref{alg. SmoothingEngineEstimation}}} }){Smoothing Engine}{
			$(\vect{\bar{X}}^{0:k+L},\vect{\Sigma}^{0:k+L}) = SmoothingEngine(\vect{\bar{X}}^{0}, \vect{\Sigma}^{0},\vect{U}^{0:k+L-1},\vect{Z}^{0:k+L}_{ML})$
			计算先验估计$b(\vect{X}^{k+1:k+L|k}) \leftarrow$ 式~\eqref{eq. SmoothingPriorEstimates}\;
		}
		\ForEach(\tcc*[f]{重复调用算法~\ref{alg. PROB_APSE}}){$s_{ij}^{k+l} \in \mathcal{E}^{k+1:k+L}_{iter}$}{
			提取$b(\vect{p}_i^{k+l}),b(\vect{X}_j^{k+l}) \leftarrow$ $b(\vect{{X}}^{k+1:k+L|k})$ \;
			$p(s_{ij}^{k+l}) = ComputeSingleConnProb(b(\vect{p}_i^{k+l}),b(\vect{X}_j^{k+l}),\rho,\delta)$
		}
		$p(\mathcal{E}^{k+1:k+L}_{iter}|\mathcal{H}^{0:k+L|k}) = \prod_{\forall s_{ij}^{k+l} \in \mathcal{E}^{k+1:k+L}_{iter}} p(s_{ij}^{k+l})  $ \;
		根据$\vect{\bar{X}}^{k+1:k+L},\vect{\Sigma}^{k+1:k+L}$的值计算目标函数$J_{\mathcal{E}^{k+1:k+_L}_{iter}} \leftarrow$ 式~\eqref{eq. ACTCO_objective_function}  \;
		\KwRet{$J_k^{iter} = J_{\mathcal{E}^{k+1:k+_L}_{iter}} \times p(\mathcal{E}^{k+1:k+L}_{iter}|\mathcal{H}^{0:k+L|k})$}\;
	}	
	\AlgoBiCaption{给定网络拓扑和候选控制时的目标函数计算算法}{Compute objective function given control candidate and network topology}
	\label{alg. EvaluateObjectiveFucntionProbCon}
\end{algorithm}

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwIn{$\rho$传感器最大测量范围； $\delta$ 连通概率计算误差阈值 \\
		\quad $b(\vect{X}^{k}) \sim \mathcal{N}(\vect{\bar{X}}^{k}, \vect{\Sigma}^{k})$:  多智能体系统在规划时刻$k$的信念分布$b(\vect{X}^{k})$\\ 
		\quad $b(\vect{X}^{0}) \sim \mathcal{N}(\vect{\bar{X}}^{0}, \vect{\Sigma}^{0})$: 多智能体的初始信念 $b(\vect{X}^{0})$  \\
		\quad $\mathcal{H}^{0:k+L|k} = \{ \mathcal{H}^{0:k}, \vect{U}^{k:k+L-1} \}$ : 包含待评估控制输入的历史数据集合 
	}
	\KwOut{\quad $J_k$: 给定$\vect{U}^{k:k+L-1}$的目标函数值 }
	
	\SetKwProg{Fn}{Function}{:}{end}	
	\Fn{ EvaluateObjective( \textbf{Input} ) } {
		\ForEach(\tcc*[f]{重复调用算法~\ref{alg. EvaluateObjectiveFucntionProbCon}}){$\mathcal{E}^{k+1:k+L}_{ter} \in \vect{E}$}{
			$J_k^{iter} = EvaluateObjectveProbCon( \textbf{Input},\mathcal{E}^{k+1:k+L}_{ter} )$ \;
		}
		$J_k = \sum_{iter} J_k^{iter}$ \;
	}	
	\AlgoBiCaption{给定候选控制时基于概率拓扑的目标函数计算算法}{Compute objective function given control candidate based on probabilistic network topology}
	\label{alg. EvaluateObjectiveFucntionAll}
\end{algorithm}


至此，对概率拓扑的信念空间规划问题，其目标函数计算过程可以从式~\eqref{eq. ACTCO_objective function1}和式~\eqref{eq. ACTCO_objective function2}总结为如图~\ref{fig. PROB_ComputeObjectiveProbCon}和算法~\ref{alg. EvaluateObjectiveFucntionAll}所示的计算框架。为了评估一组候选的控制输入$\vect{U}^{k:k+L-1}$，需要计算全部可能拓扑分支对应的目标函数值，然后将拓扑概率作为权值，对目标函数值进行权值累加。
\begin{figure}[ht]
	\centering
	\includegraphics[width = 0.95\textwidth]{figures_PROB/computeObjective.jpg}
	\bicaption[fig. PROB_ComputeObjectiveProbCon]{}{基于概率拓扑的目标函数计算示意图}{Fig.$\!$}{The computation of objective function under probabilistic network}
\end{figure}

算法~\ref{alg. EvaluateObjectiveFucntionAll}是对前一章目标函数计算方法~\ref{alg. EvaluateObjectiveFucntion}的改进。因此，基于概率拓扑的构型规划方法是在前一章的规划框架( 算法~\ref{alg. EvaluateObjectiveFucntion}-\ref{alg. PlanningScheme} )基础上，使用基于概率拓扑的目标函数计算算法~\ref{alg. EvaluateObjectiveFucntionAll}代替算法~\ref{alg. EvaluateObjectiveFucntion}，其他部分均保持不变。


\BiSection{数值仿真结果}{Simulation results}

本节将通过两组数值仿真实验验证本章所提方法的正确性和有效性，分别为单智能体节点观测连通维持问题(Connection Maintenance, CM)以及多智能体系统中基于主动规划的协同定位问题(Cooperative Localization using Active Planning, CLAP)。具体的场景配置将在后文给出。

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwOut{ $p(s_{1,2}^{k+1}=1)$  }
	\KwIn{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}},\rho$ \;}
	
	
	\SetKwProg{Fn}{Function}{:}{end}	
	\Fn{ BernoulliModel( \textbf{Input} ) } {
		$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
		$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
		
		$p(s_{1,2}^{k+1}=1) = 0$ , $d = ||\vect{\mu}_{\Delta}||_2$\;
		\If{$d<\rho$}{
			$p(s_{1,2}^{k+1}=1) = 1$ \;
		}
		
		\textbf{Return} $p(s_{1,2}^{k+1}=1) $\;
	}
	\AlgoBiCaption{伯努利概率模型}{The Bernoulli model}
	\label{alg. BernoulliDistributionAlg}
\end{algorithm}
\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwOut{ $p(s_{1,2}^{k+1}=1)$  }
	\KwIn{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}},\rho$ \;}
	
	
	\SetKwProg{Fn}{Function}{:}{end}	
	\Fn{ LinerModel( \textbf{Input} ) } {
		$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
		$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
		
		$[ \rho_{3\delta}^{min}, \rho_{3\delta}^{max} ] \leftarrow \Delta \vect{p}_{12}^{k+1}$节点分布的$3\delta$区域计算(参见\text{图~\ref{fig. PROB_3delta} }) \;
		
		\uIf{$\rho > \rho_{3\delta}^{max}$}{
			\textbf{Return} $p(s_{1,2}^{k+1}=1) = 1$\;
		}
		\ElseIf{$\rho < \rho_{3\delta}^{min}$}{
			\textbf{Return} $p(s_{1,2}^{k+1}=1) = 0$\;
		}
		$k = \frac{1}{\rho_{3\delta}^{max} - \rho_{3\delta}^{min}}$ \;
		
		\textbf{Return} $p(s_{1,2}^{k+1}=1) = k(\rho - \rho_{3\delta}^{min})$\;
	}
	\AlgoBiCaption{线性概率模型}{The Linear model}
	\label{alg. LinearDistributionAlg}
\end{algorithm}

本章的主要创新点在于，设计了一类圆盘观测模型连通概率的预测算法APSE，通过严格的理论推导给出了概率预测误差小于给定阈值的充分条件，并将该算法应用于构型规划问题的目标函数评估。已有文献对圆盘观测模型的概率预测方法多是基于经验的近似方法，本章将把APSE算法作为实验组，与三种近似概率计算方法进行对比，分别为：
\begin{enumerate}
	\item[(1)] \textbf{伯努利模型：}即上一章中基于最大似然观测假设的概率计算方法。未来时刻的连通变量是一类伯努利变量，其概率值仅可取0或1。具体计算过程如算法~\ref{alg. BernoulliDistributionAlg}，该算法仅使用两个节点相对分布的均值向量与观测传感器最大覆盖范围进行比较，从而确定两个节点之间的观测连通存在与否。
	\item[(2)] \textbf{线性模型：}文献\cite{indelman2015planning}使用的一种连通概率近似计算方法，观测连通存在与否的概率与节点之间相对距离成反比例关系。根据圆盘模型的分布特征，本节对文献\cite{indelman2015planning}中的线性模型进行了改进，加入了$3\delta$的校正方法，具体计算如算法~\ref{alg. LinearDistributionAlg}所示。
	\item[(3)] \textbf{随机抽样模型：}即小节~\ref{subsec. ch5sec5.3.5}中使用的对比算法~\ref{alg. RandomSampleDistributionAlg}。由于前文的仿真对照实验已经表明，本章所设计的APSE算法的计算复杂度位于抽样规模为$dim=10^3-10^4$的随机抽样算法之间。因此作为对照组，本节设定随机抽样算法的抽样规模为$dim=10^4$。
\end{enumerate}

对照组和实验组在构型规划框架中的唯一区别位于算法~\ref{alg. EvaluateObjectiveFucntionProbCon}第$16$行，将分别调用四种算法对给定的未来观测拓扑分支$\mathcal{E}^{k+1:k+L}_{iter}$中每一个观测连通的概率进行计算。

另外，当规划框架采用随机模型计算连通概率时，需要对所有可能拓扑的概率进行归一化处理。APSE、伯努利模型以及线性模型都是确定性算法，即当算法输入不变，多次运行的输出也不变。然而随机抽样模型是不确定性算法，即使算法的输入是一样的，算法~\ref{alg. RandomSampleDistributionAlg}运行两次得到的输出会因为随机数生成过程的不确定性而稍稍有所差别。由于所有可能拓扑集合$\vect{E}$中两个相邻元素通常有许多相同的观测连接，即$\mathcal{E}^{k+1:k+L}_{iter} \bigcap \mathcal{E}^{k+1:k+L}_{iter+1} \neq \varnothing$，而在算法~\ref{alg. EvaluateObjectiveFucntionProbCon}中这些相同观测连接的概率需要被重复计算，因此有可能导致$\sum_{\forall iter} p(\mathcal{E}^{k+1:k+L}_{iter}) \neq 1$。因此当使用随机抽样模型时，在计算得到所有可能拓扑的概率$p(\mathcal{E}^{k+1:k+L}_{iter})，\forall iter$后，需要对它们进行归一化处理，然后才能作为权值对$J_{iter}^k$进行累加。

\BiSubsection{CM仿真案例}{The simulation case of CM}

本小节将详解介绍并分析CM仿真的配置与结果。本组仿真假定在GNSS拒止的环境下存在一个静止的锚节点$A$以及一个可以移动的智能体$B$。锚节点知道自身的位置信息，并装备有通信、观测传感器，可以对传感器覆盖范围$\rho$内的节点进行测量，并持续广播自身位置。智能体仅知道其初始位置，并且也配备了与锚节点相同的通信、观测传感器。智能体的运动和观测过程均存在噪声，因此智能体在每次运动后都需要将内部传感器数据与对外部环境的观测数据进行融合从而进行自定位。假初始时刻智能体位于锚节点观测通信区域外围。由于环境中不存在其他的可观测信息，因此智能体的运动目标是首先与锚节点获取观测连接，然后根据与锚节点的观测校正自定位误差，并通过规划自身运动轨迹在最小化自定位误差的前提下，与锚节点之间保持最大的相对距离。该任务配置是GNSS拒止环境下的一类典型巡视任务场景，在缺乏定位信息的环境中，智能体需要在锚节点辅助的情况下，在自定位的同时，尽可能地对锚节点未能覆盖的区域进行巡视覆盖。

使用$\vect{\bar{p}}_{B}^{k+l}$和$\vect{\Sigma}_{\vect{\bar{p}}_{B}^{k+l}}$表示智能体在第$l$个规划时刻对自身位置状态的预测均值和协方差。则式~\eqref{eq. ACTCO_objective_function}中的瞬时目标函数为：
$$c^l(\cdot) = \omega_1 tr(\vect{\Sigma}_{\vect{\bar{p}}_{B}^{k+l}}) + \omega_2 |\bar{d}_{AB} - \rho| $$
式中$\bar{d}_{AB} = || \vect{\bar{p}}_{B}^{k+l} - \vect{\bar{p}}_{A}^{k+l} ||_2$为智能体对其与锚节点之间相对距离的预测。$\omega_1,\omega_2$为两部分目标的权值。从目标函数的组成上可以看出，智能体$B$的最优策略是一直保持在锚节点观测范围的边界运动。

\begin{figure}[htbp]
	\centering	
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Bernoulli model]{\subfigure[伯努利模型]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/1CMDistributionOver50Steps1}}}	
	\null\hfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Linear model]{\subfigure[线性模型]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/2CMDistributionOver50Steps1}}}
	\vfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Random Sampling model]{\subfigure[随机抽样模型]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/3CMDistributionOver50Steps1}}}	
	\null\hfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[APSE]{\subfigure[APSE]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/4CMDistributionOver50Steps1}}}	
	
	\bicaption[fig. PROB_performanceRepresentativeCase]{}{CM案例四种方法的性能对比}{Fig.$\!$}{Performance comparison for the CM scenario}
\end{figure}
\begin{figure}[htbp]
	\centering	
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Mean of leader's localization uncertainty over all trials]{\subfigure[领导者定位的平均不确定性]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/4CMMeanFourInOneOver50Steps2}}}	
	\null\hfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Mean of leader's localization uncertainty over all trials and over all time steps as well as the totall connections]{\subfigure[领导者定位双层平均误差以及总连接个数]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/NumConAndMeanUncertainty}}}
		
	\bicaption[fig. PROB_performanceRepresentativeCase]{}{CM仿真四种方法的数值结果对}{Fig.$\!$}{Statistical results and comparsion for the CM scenario}
	\vspace{-0em}
\end{figure}

本节使用前一章所设计的基于数值梯度的信念空间构型规划框架，并将目标函数计算模块替换为本章的算法~\ref{alg. EvaluateObjectiveFucntionProbCon}-\ref{alg. EvaluateObjectiveFucntionAll}。信念空间构型规划框架的参数与前一章仿真中的参数保持相同，参见表~\ref{tab. ParametersFreeSpace}。锚节点和智能体的初始位置为$\vect{P}_A^0 = [0,0]^T,\vect{P}_B^0 = [8,2]^T$。节点配备的通信、观测传感器最大测量范围为$\rho= 5 m$，因此在初始时刻$||\vect{P}_A^0 = [0,0]^T-\vect{P}_B^0||_2 > \rho$，即初始时刻两个节点之间不存在观测连通。算法~\ref{alg. PROB_APSE}中连通概率的期望误差阈值$\delta$设置为$10^{-10}$，规划期限$L=5 s$，目标函数的增益权值为$\omega_1=0.9$,$\omega_2=0.1$。三个对照组方法和APSE实验组的方法分别重复仿真30次，每次持续50步。选用智能体节点的自定位误差为衡量指标，并采用节点自定位分布协方差的迹作为自定位误差。在30次打靶仿真中，智能体使用四种概率计算方法得到的自定位误差随时间变化的分布如图~\ref{fig. PROB_performanceRepresentativeCase}所示，四个子图的纵坐标刻度是相同的。

针对测试的四组方法，智能体自定位误差在30次打靶仿真里随时间增长平均变化曲线的对比图可参见图~\ref{fig. PROB_performanceRepresentativeCase} a)。在总共1500个仿真时刻中( 30次打靶仿真，每次打靶仿真持续50个时刻 )，智能体节点的自定位误差均值，以及智能体节点与锚节点之间总共存在连接次数的统计如图~\ref{fig. PROB_performanceRepresentativeCase} b)所示。


首先，从仿真结果的对比图中可以看出，四组概率预测方法都得出了相同变化趋势，即智能体自定位误差随时间的增长而增长。这是因为，智能体无法保证保证一直处于锚节点的观测范围之内，此时它的自定位过程将不存在外部绝对信息的校正，然而运动噪声却始终在累加。因此智能体节点自定位误差的长期发散趋势是正确的。

其次，自定位误差的增长速度按照伯努利模型、线性模型、随机抽样模型以及APSE方法的顺序依次减少，恰好对应于四种方法在预测连通概率上预测精度由低到高的顺序。因此可以总结，在规划期限内观测拓扑的概率预测精度越高，节点自定位误差增长约缓慢，即定位效果越好。

最后，图~\ref{fig. PROB_performanceRepresentativeCase} b)的数据可以解释概率预测精度对协同定位的影响，即未来拓扑概率预测地约精准，智能体越倾向于与锚节点建立更多的观测连接。该结论与已有成果中观测越多定位效果越好的结论是一致的。

\BiSubsection{CLAP仿真场景}{The simulation case for CLAP}

本小节将详解介绍和分析CLAP任务的仿真配置与结果。如图~\ref{fig. PROB_CLAPComposition}所示，在GNSS拒止的自由任务空间中，存在一个锚节点和两个智能体节点，分别为领导者$L$以及跟随者$F$。其中跟随者的主要任务是通过规划自身路径为领导者提供中继的观测和通信连接，从而帮减少领导者的自定位误差。领导者是一类具有预定义高层任务的节点，由于它的预设任务轨迹完全位于锚节点的观测通信范围之外，因此需要使用跟随者节点作为中继，为领导者提供外部观测数据以校正其自定位误差。当跟随者节点同时与领导者，以及锚节点建立了观测连接后，在某种意义上等价于领导者与锚节点之间存在了虚拟观测连接。该虚拟观测连接对领导者定位误差校正作用的效果取决于跟随者、领导者以及锚节点所组成多智能体系统的构型。又由于这三个节点中只有跟随者节点的位置是完全可控的，因此CLAP问题的核心是规划跟随者节点的实时位置。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/sim/CLAPinitialPlot}
	\bicaption[fig. PROB_CLAPComposition]{}{CLAP问题的场景配置}{Fig.$\!$}{The composition for CLAP problem}\vspace{0em}
\end{figure}

本节使用$\vect{\bar{p}}_F^{k+l},\vect{\Sigma}_{\vect{\bar{p}}_F^{k+l}}，\vect{\bar{p}}_L^{k+l},\vect{\Sigma}_{\vect{\bar{p}}_L^{k+l}}$表示跟随者节点在第$l$个规划时刻对其自身位置状态和领导者位置状态的预测均值和协方差矩阵。则式~\eqref{eq. ACTCO_objective_function}中的瞬时目标函数为：
$$c^l(\cdot) = \omega_1 tr(\vect{\Sigma}_{\vect{\bar{p}}_L^{k+l}}) + \omega_2 tr(\vect{\Sigma}_{\vect{\bar{p}}_F^{k+l}})$$

领导者和跟随者的初始位置分别为$\vect{{p}}_F^{0} = [6,0]^T$以及$\vect{p}_L^0 = [3,2]^T$并且有一定的定位误差，即图~\ref{fig. PROB_CLAPComposition}中领导者和跟随者节点周围的误差椭圆。锚节点被放置于原点。初始时刻，跟随者位于锚节点的观测通信范围内，并且与领导者之间也存在观测连接。领导者给定的高层任务轨迹由四个任务点所组成，分别为$W_1 = [6,6]^T,W_2 = [-6,6]^T,W_3 = [-6,-6]^T,W_4 = [6,-6]^T$。领导者节点每一时刻都根据对自身位置状态的估计$\vect{\bar{p}}_L^k$以及任务点$W_i$，按照下式计算其控制输入，
$$\vect{u}_L = u_{max} \frac{\vect{\bar{p}}_L^k - W_i}{||\vect{\bar{p}}_L^k - W_i||_2}$$
式中$u_{max}$为智能体节点每个时间步长能行走的最大路程，本小节的CLAP仿真设置领导者的$u_{max}=0.5 m$，跟随者的$u_{max}=0.7 m$。并且当$||\vect{\bar{p}}_L^k - W_i||_2<0.2 m$时，领导者将依次切换它的目标任务点。瞬时目标函数中的增益为$\omega_1 = 9, \omega_2=1$。为了验证使用各种概率计算算法的构型规划方法对环境不确定性的鲁棒性，CLAP问题对四种概率计算方法各自进行20次打靶仿真，以总结表现性能的概率特征。

在20次打靶中，领导者节点对应于四种概率预测算法的自定位误差随时间变化的分布如图~\ref{fig. ACOTO_performanceCLAP4}所示，其中四张子图的纵坐标尺度是\textbf{不一致}的。
与CM问题不同，CLAP问题的结束是当领导者从初始位置出发，沿着给定的任务点序列运动，直到到达最后一个任务点时整个任务结束。由于受到随机噪声的干扰，即使在使用同一种概率预测方法的两次打靶仿真中，领导者在同一时刻的自定位也会有所不同。因此在所有四组方法的20次打靶仿真中，每一次打靶仿真持续的仿真时间都不一样。
因此图~\ref{fig. ACOTO_performanceCLAP4}在每个子图的垂直方向上增加了四个阴影区，分别对应领导者节点到达每一个任务点的时间分布。每个子图顶端的颜色图表征了每组方法的20次打靶仿真在横轴对应的时间上，还有多少次尚未结束(即领导者节点尚未到达最后一个任务点)，右端的颜色图给出了尚未结束仿真的数目。

在20次打靶仿真中，四组方法得到的领导者节点自定位误差的均值分布如图~\ref{fig. ACOTO_performanceCLAP2} a)所示。领导者节点的估计位置与真实位置之间误差的RSME平均分布如图图~\ref{fig. ACOTO_performanceCLAP2} b)所示。


表~\ref{tab. performanceCLAP}是四组实验20次打靶仿真中四个通用指标的统计情况。这四个指标分别为(指标名称中第一个“m”表示取20次打靶仿真里该指标的平均值，第二个“m”表示取每一个仿真中所有仿真时刻里该指标的平均值。)

\begin{figure}[htbp]
	\centering	
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Bernoulli model]{\subfigure[伯努利模型]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/CLAPUncertaintyDistribution1_20}}}	
	\null\hfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Linear model]{\subfigure[线性模型]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/CLAPUncertaintyDistribution2_20}}}
	\vfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Random Sampling model]{\subfigure[随机抽样模型]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/CLAPUncertaintyDistribution3_20}}}	
	\null\hfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[APSE]{\subfigure[APSE方法]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/CLAPUncertaintyDistribution4_20}}}
	
	\bicaption[fig. ACOTO_performanceCLAP4]{}{CLAP案例四种方法的性能对比}{Fig.$\!$}{Performance comparison for the CLAP scenario}
	\vspace{-0em}
\end{figure}

\begin{enumerate}
	\item[(1):] \textbf{mm-RMSE:} 领导者节点估计位置与真实位置之间误差RMSE的双层平均值(对20次打靶仿真和每一次打靶仿真中的所有时刻里的该指标求取平均值)。
	\item[(2)] \textbf{mm-Trace:} 领导者节点自定位误差的双层平均值，采用求自定位协方差迹来表示自定位误差。
	\item[(3)] \textbf{TotalCon:} 每一种方法在所有20次打靶仿真以及所有仿真时刻中，领导者与锚节点存在虚拟连接的总数。
	\item[(4)] \textbf{m-TotalTra:} 领导者节点从初始位置到最后一个任务点所行走总路程的单层平均值(即将20次打靶仿真得到的该指标求取平均值)。
\end{enumerate}

\begin{figure}[htbp]
	\centering	
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Mean RMSE for the deviations of leader's estimated and real position]{\subfigure[估计位置与真实位置的平均偏差RMSE]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/CLAPRMSE2}}}	
	\null\hfill
	\subfigure{}\addtocounter{subfigure}{-2}
	\subfigure[Mean of leader's localization uncertainty over all trials]{\subfigure[领导者的平均自定位误差]{\includegraphics[width = 0.48\textwidth]{figures_PROB/sim/CLAMPMeanUncertaintyDistribution}}}
	
	\bicaption[fig. ACOTO_performanceCLAP2]{}{CLAP案例四种方法数值结果对比}{Fig.$\!$}{Statistical performance comparison for the CLAP scenario}
	\vspace{-0em}
\end{figure}


从这些对比结果中可以总结处两点结论：一，基于APSE方法的构型规划方法在减少领导者自定位误差方面的表现优于其他所有对照组；二，构型规划中使用的连通概率算法的预测精度越高，领导者自定位误差越小，在20次打靶仿真中的分布越集中，随时间增长的速度越慢。经计算发现，采用本章设计的APSE算法进行构型规划实验组，比使用\{ 伯努利模型，线性模型，随机抽样模型 \}的方法的对照组，在mm-RMSE指标上提升了\{ 24.28\%,22.18\%,11.40 \} 在mm-Trace指标上提升了\{ 56.86\%, 42.36\%, 17.74\% \}；在TotalCon指标上提升了\{ 81.24\%, 75.73\%, 13.47\% \}。

\begin{table}[htbp]
	\bicaption[tab. performanceCLAP]{}{CLAP问题的四个指标统计情况}{Table$\!$}{Statistical characteristics summarized over 20 trials and over all simulating time steps.}
	\centering
	\wuhao		
	\begin{tabular}{ccccccc}
		\toprule[1.5pt]
		\multirow{2}{*}{性能指标} & \multicolumn{4}{c}{ $BSP$规划框架 } \\ \cline{2-5}
		& 伯努利模型 & 线性模型 & 随机抽样概率($dim=10^{4}$) & APSE \\ 
		\midrule[1pt]
		mm-RMSE (m)              & 1.0669   & 1.0281   & 0.9119    & \textbf{0.8079}    \\
		mm-Trace (m)            & 1.5651   & 1.1715   & 0.8208    & \textbf{0.6752}    \\
		TotalCon              & 837  & 894  & 1337   & \textbf{1517}   \\
		m-TotalTra (m)             & 47.1090  & 46.9427  & 46.6451   & $\mathbf{45.7925}$   \\ 
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

对于最后一项指标，如果将该指标与领导者节点给定路径的标称路程进行比较，则可以反应运动过程和测量过程的噪音在领导者位置状态上的累计效应。通过计算，领导者节点的标称路程是$42 m$，四组方法(按照表~\ref{tab. performanceCLAP}排序)因扰动而累计的额外路程是\{ 5.1090, 4.9427, 4.6451, 3.7925 \}m. 因此实验组在该项指标上比三个对照组的提升分别为\{ 25.76\%, 23.27\%, 18.35\% \}。



\BiSection{小结}{Summary}

本章以前一章设计的信念空间构型规划方法为基础，进一步研究规划期限内未来观测不可知问题对目标函数计算过程的影响，提出了一类基于拓扑概率的构型规划方法。该方法与前一章方法的最大不同之处在于，给定一组候选控制后，前一章的方法仅使用有最大似然观测的未来分支来评估该控制输入的性能；然而基于概率拓扑的规划方法将评估每一个可能观测拓扑的未来分支，并将拓扑概率作为权值，通过加权计算来评估该控制输入的性能。本章首先将一类圆盘观测模型的连通建模为一类正态二次型变量，并使用有限项近似的多项式序列展开方法计算圆盘观测模型的连通概率。通过严格的理论推导，本章给出了截断误差稳定的充分条件，从而从理论上保证了圆盘模型未来观测连通概率的预测精度，总结了圆盘观测模型的概率预测算法APSE。然后本章将APSE算法应用到信念空间的构型规划问题中，提出了基于概率拓扑的构型规划框架。最后通过两个场景下的构型规划任务仿真，验证了本章所提方法的正确性和有效性，得出了如下结论：

\begin{enumerate}
	\item[(1)] 本章针对圆盘观测模型所设计的自适应概率计算方法APSE具有理论保障的计算精度，单次运算时间小于1.5ms，计算复杂度介于$10^{3}-10^{4}$的抽样算法之间，然而概率计算精度与$10^{7}$的抽样算法相同。
	\item[(2)] 基于概率拓扑的构型预测算法对多智能体系统的未来预测更为精准，得到的规划方案对环境不确定性具有更强的鲁棒性，且拓扑概率越精确鲁棒性越强。
\end{enumerate}

本章所提方法存在的问题主要有：
\begin{enumerate}
	\item[(1)] APSE算法需要对通信距离$\rho$进行大量的幂次计算，因此当$\rho$的取值较大时，计算平台的舍入误差将对APSE算法的计算精度产生较大影响。
	\item[(2)] 基于概率拓扑的构型规划方法具有指数增长的计算复杂度。本章设计的方法需要评估所有可能的观测拓扑，而未来观测拓扑的数量随着多智能体系统规模以及规划期限长度的增加而指数增长。该问题有待进一步研究。
\end{enumerate}
